Atomistic simulation of structure and dynamics of columnar phases of 

hexabenzocoronene derivatives 
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Using atomistic molecular dynamics simulations we study solid and liquid crystalline columnar 
discotic phases formed by alkyl-substituted hexabenzocoronene mesogens. Correlations between the 
molecular structure, packing, and dynamical properties of these materials are established. 



I. INTRODUCTION 

Most discotic thermotropic liquid crystals are formed 
by flat molecules with a central aromatic core and sev- 
eral aliphatic chains attached at the edges 0, 0. The 
size and shape of the core can be varied, as well as the 
length and structure of the side chains, which helps to 
control the functional properties of these materials Q. 
Triphenylenes Jj] , hexabenzocoronenes (HBC) and their 
derivatives d 0] are the most frequently studied exam- 
ples of discotics. 

The main feature of discotic liquid crystals is that 
they form columnar phases, where the molecules stack 
on top of each other and the columns arrange in a regu- 
lar lattice as sketched in Fig.^ The self-organization 
into stacks with aromatic cores surrounded by saturated 
hydrocarbons results in predominantly one-dimensional 
charge transport within the core, along the columns 0,0- 
By modifying the periphery and the shape of the core it 
is possible to obtain materials with small charge trapping 
and recombination. 

Conducting properties of discotics are useful for var- 
ious photovoltaic devices, in particular solar cells and 
plastic organic field-effect transistors. In solar cells it is 
important that the charges can find a direct percolation 
path to the electrodes without being trapped in an is- 
land of the material. Therefore, the morphology of the 
conducting film has a big impact on the efficiency of the 




FIG. 1: Columnar phase formed by the disc-shaped molecules 
and the most common arrangements of columns in two- 
dimensional lattices: (a) hexagonal, (b) rectangular, and (c) 
herringbone. 



device [IHEJ. The ideal structure is an aligned discotic 
liquid crystal with the columns perpendicular (solar cell) 
or parallel (field-effect transistor) to the substrate. 

The spacial arrangement of stacks is, however, never 
perfect: the columns can be misaligned, tilted, or form 
various types of topological defects. In addition, the local 
alignment of molecules in columns can be different for dif- 
ferent compounds, changing the overlap of the 7r orbitals 
and affecting the efficiency of the charge transport in a 
single column. It is, therefore, important to control both 
the local alignment of the molecules in the columns and 
the global arrangement of the columns in the mesophase. 
Since these properties are fixed by the molecular architec- 
ture and processing conditions, the main task of synthetic 
chemistry is to synthesize compounds which have various 
shapes of the core and different peripheries, to optimize 
the conducting properties of the mesophase. One should 
keep in mind, however, that, perfect arrangement of the 
molecules itself does not always mean optimal conduct- 
ing properties: overlap of the ir — ir orbitals, splitting of 
the HOMO level, etc. affect the efficiency of the charge 
transport and should be taken into account. 

The actual synthesis of each particular compound can 
involve many stages and be rather time-consuming. It is 
thus helpful to get an insight into the mesophase proper- 
ties prior to the actual synthesis. For this purpose various 
computer simulation techniques can be used. Depending 
on the length- and time-scales several techniques can be 
involved: for example quantum chemistry calculations 
are used to study the adsorption of molecules on elec- 
trodes and the charge transfer mechanisms ; MD sim- 
ulations with atomistic force fields are employed to calcu- 
late local properties, such as order parameters and molec- 
ular arrangement [lj, [FJLLJ UH j finally, coarse-grained 
simulations are successful in studying the morphology of 
the bulk material, global arrangement of columns, generic 
phase diagrams, and defects [H E 0, . 

The number of computer simulations so far performed 
on discotic materials, is, however, relatively small and 
most of them are done using idealized model poten- 
tials 0, 0]. However, a closer contact with real ex- 
periments requires more chemical details to be put in. 
Up to now, there have only been a few attempts to per- 
form atomistic simulations of bulk discotic liquid crystals, 
mostly on the triphenylene-based systems with a small 
aromatic core I5j. Derivatives of hexabenzocoronenes 
are, however, more promising for molecular applications, 
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FIG. 2: Studied derivatives of hexabenzocoronene: R = Cio, 
Ci 2 , C14, Ci 6 , Cio-6, PhCi 2 . Only R = C12 , C10-6 and 
PhCi2 are shown. 



due to their larger aromatic core and better overlap of 
the 7r-orbitals. The atomistic simulation of the structure 
and dynamics of the derivatives of hexabenzocoronenes 
is the subject of the present paper. 



II. STUDIED SYSTEMS AND 
COMPUTATIONAL DETAILS 

The molecular structures of the studied hexabenzo- 
coronene derivatives are shown in Fig. [21 All of them con- 
sist of a central flat aromatic core and six side chains. We 
considered several types of side chains: alkyl chains of dif- 
ferent lengths, C n , with n — 10, 12, 14, 16; branched side 
chains, Cio-e! an d dodecylphenyl-substituted PhCi2. 
The choice of systems is based on the availability of ex- 
perimental data (wide angle X-ray scattering, solid-state 
NMR, and differential scanning calorimetry 0, I2J 0, 

mil). 



A. Model potential 




FIG. 3: Molecular fragment which illustrates the atom types 
used in the MD simulations. Methylene and methyl groups 
of the aliphatic chains are considered as one interaction site, 
according to the united atom approach. 



In our simulation we adopt the united atom approach, 
and consider explicitly only the hydrogen atoms belong- 
ing to the aromatic rings of the central core, since the 
charge on these hydrogens is responsible for the multi- 
pole of the core of the molecule and the interaction of 
the side chains with aromatic cores is not important in 
the columnar phase. To describe the inter- and intra- 
molecular interactions, we define eight different types of 
atoms (see Fig. |3] and Table QJ: aromatic carbons of the 
central ring (Cb), aromatic carbons bonded to the hydro- 
gens (Ch), hydrogens (H), aliphatic carbons representing 
either a methylene group (C2), terminal methyl group 
(C3) or a junction (Ci), and aromatic carbons bonding 
either alkyl chain (Ca) or two aromatic subsystems (Cp). 

The potential energy of the system is then given by the 
following sum of contributions 
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The first four terms in Eq. ^ correspond to the energy 
due to the bonded interactions: harmonic bonds (1-2 in- 
teractions), angles (1-3 interactions), and torsions, which 
include Ryckaert-Bellemans and improper dihedrals (1-4 
interactions). The parameters of the bonded potentials 
are taken from the AMBER ^ a nd OP LS 0, HJ 
force fields and are listed in Tables ITllIVI 

The last double sum in Eq. ^is due to the nonbonded 
interactions. It is given by the Coulomb and Lennard- 
Jones terms with qi the effective charges and and 
Uij the well depth and contact distance of the sites i 
and j. The interaction potential between two atoms of 
the same type has its own set of parameters {q, e, a}, 
whereas between two different atoms the interaction en- 
ergy is given by the geometrical mean of the two pairs 
of atoms and the arithmetic average for the lengths a 
(Lorentz-Berthelot mixing rules). Note that 1-2, 1-3, 
and also 1-4 interactions do not contribute to the above 
double sum. Furthermore, nonbonded interactions be- 
tween pairs of sites belonging to the central aromatic 
core have been excluded. The parameters of the non- 
bonded interactions are given in Table |T| Note that 
a similar model has already been used to study hex- 
akis(pentyloxy)triphenylenes |15| . 
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atom 


mass (au) 


9(e) 


e(kJ-mor 1 ) 


<r(nm) 


Ca, Cb, Cp 


12.011 





0.293 


0.355 


Ch 


12.011 


-0.115 


0.293 


0.355 


H 


1.008 


0.115 


0.125 


0.242 


Ci 


13.019 





0.450 


0.3905 


c 2 


14.027 





0.494 


0.3905 


c 3 


15.035 





0.732 


0.3905 



TABLE I: Atom types and parameters of nonbonded interac- 
tions taken from the OPLS force field. Atom types correspond 
to those shown in Fig. [3 



bond 


b (nm) 


K^kJ-moP 1 - nm- 2 ) 


C[A,B,H,P] — C[A,B,H,P] 


0.139 


418400 


C[2,3] — C[ 2j 3] 


0.153 


334720 


Ch-H 


0.108 


418400 


Ci — C2 


0.153 


334720 



TABLE II: Force constants and equilibrium bond distances. 
Atom types correspond to those shown in Fig. |3] Parenthe- 
ses mean that no differences have been made between the 
included types of atoms. 



angle 



C[B,H] — C[b,H,2] — C[B,H,2] 
C2 — C2 — C[2,3] 

Cb — C2 — Ci 
Ci — C2 — C2 



g(deg) 



120 
109.47 

120 
109.47 



Kg (kj ■ mol • rad 



527 
502.1 

527 
502.1 



TABLE III: Force constants and equilibrium angles. Atom 
types correspond to those shown in Fig. Parentheses mean 
that no differences have been made between the included 
types of atoms. 



B. Computational details 

We simulated systems of 160 molecules, arranged in 
columns of 10 molecules each. The initial configura- 
tion was either a hexagonal (T = 400 K) or rectan- 
gular (T — 300 K) arrangement of the columns (see 
FigCJ. The distance between the columns and between 
the molecules in the columns has been slightly increased 
to avoid site superposition. After the energy minimiza- 
tion, a short equilibrating run (4 ns) was performed at 
a constant pressure P = 0.1 MPa and a constant tem- 
perature T — 200 K. Quenching the system at a lower 
temperature helps to avoid trapping it in a disordered 
(non-equilibrium) state from the very beginning. We 
used Berendsen thermostat with anisotropic pres- 
sure coupling. The simulation box angles were fixed at 
90deg. The electrostatic interactions were treated with 
the smooth particle-mesh Ewald method [3^ (spacing for 
PME grid 0.12 nm, cutoff 0.9 nm). 

After the equilibration, the simulation box was a par- 
allelepiped of edges L x , L y , and L z , with L x /L y = y/3/2 
or L x /L y = 1, depending on the arrangement of the 
columns (hexagonal or rectangular lattice). 

The production run (100 ns) was performed at a con- 
stant pressure (P = 0.1 MPa) and temperature fixed us- 
ing Berendsen thermostat with time constant r p = 5ps 
and compressibility 4.5 x 10 _5 Bar~ 1 (both isotropic). 
Two representative temperatures were chosen: T = 
400 K and T — 300 K. The time step of 2 fs was used to 
integrate equations of motion using the velocity Verlet 
algorithm with a 0.9 nm cut-off for short-range interac- 
tions. All calculations were performed with the parallel 
version of the GROMACS program [H|. Bond lengths 
were constrained using the LINCS algorithm Q pro- 
vided by GROMACS. 



dihedral 


Vi 


v 2 


v 3 


X — C[A,B,P,H] — C[A,B,P,H] — X 


0.0 


121.34 


0.0 


C[B,2] — C[i,2] — C2 — C[2,3] 


5.9 


-1.14 


13.16 


improper 


K d 


ipo 




Cp — Ch — Ch — Cp 


167.4 







Ch — Ch — Ca — H 


167.4 







Ch — Cp — C[h,bj — H 


167.4 







CI — C2 — C2 — C2 


334.8 


35.26 





TABLE IV: Parameters of the torsional potential. X can be 
an aromatic carbon or a carbon bonded to a hydrogen. Paren- 
theses mean that no differences have been made between the 
included types of atoms. Vi, V2, and V3 are given in kj-mol -1 
for Ryckaert-Bellemans dihedrals and in kj ■ mol -1 • rad -2 
for improper dihedrals. The main task of improper dihedrals 
with Ch atom type is to keep the corresponding atoms in a 
plane. The CI — C2 — C2 — C2 improper constrains the ter- 
tiary carbon CI to tetrahedrality. Atom types correspond to 
those shown in Fig. 



III. RESULTS AND DISCUSSION 

An equilibrated snapshot of one of the systems at 
T = 400 K is shown in Fig. 21 One can see that the aro- 
matic cores are well aligned and the columns are arranged 
in an almost perfect hexagonal lattice, with the axially 
symmetric distribution of the side chains. At T = 300 K 
the columns are arranged in a rectangular lattice and 
the cores are also well aligned. To further quantify the 
arrangement of the molecules we calculated their orienta- 
tional ordering as well as several positional distribution 
functions. 



A. Orientational order 

To calculate the orientational order parameter of the 
columnar mesophase we considered only the aromatic 
cores of the molecules. The plane of the core was speci- 
fied by the six vectors connecting the six Ca atoms. The 




FIG. 4: MD simulation results: snapshot of the system with 
the C12 side chains. Aromatic cores are highlighted. Both 
top and side views are shown. T = 400 K, P = 0.1 MPa. 

vector was then denned as an average of the six nor- 
mals obtained by the cross product of every pair of two 
neighboring Ca — Ca vectors. We calculated the orien- 
tational order tensor for each column and for the whole 
system 

where is a unit vector normal to the i-th aromatic 
core, N is the number of molecules in a single column or 
in the whole system; correspondingly, the summation is 
performed either over a single column or over the whole 
system. < • • • > denotes the time average. Diagonalizing 
this tensor we obtain the order parameter Q, which is the 
largest eigenvalue of Q a f3, and the average orientation of 
the molecules, or director n (corresponding to the eigen- 
vector Q). Q = 1 implies perfect alignment of molecules 
in the column (all unit vectors parallel to each other), 
Q = corresponds to the isotropic angular distribution 
of the unit vectors. 

The order parameter Q and the average tilt angle of the 
unit vector n (given by the projection onto the z-axis) 
are summarized in FigEl It is clear that the systems with 
linear side chains Cio — C\§ are well ordered: the order 



FIG. 5: Order parameter Q and the ^-component of the di- 
rector (average tilt in the column) as a function of column 
number. T = 400 K. 

parameter is close to its maximum value, Q = 1; in addi- 
tion, there is no tilt of the molecules in the columns. The 
system with the dodecylphenyl-substituted side chains, 
PhCi2, has slightly lower ordering (but still rather high, 
Q « 0.95) and a small average tilt of the molecules in the 
columns. Finally, the system with the branched chains, 
Cio-6, has the worst ordering: the order parameter is 
quite low and varies much from column to column; there 
is also a significant tilt of the molecules in the columns. 
As can be seen in Table fvTI in the case of linear chains, 
the values of the order parameter Q at T — 300 K are 
identical to the one at T = 400 K , while the PhCi2 shows 
a bigger disorder; the Cio-e becomes more ordered with 
the decrease in temperature. 

B. Positional order 

To characterize the bulk structure of the columnar 
phase, we have calculated two additional correlation 
functions. The first one, <7 z (r), provides information 
about the arrangement of the molecules in the planes 
parallel to the director n, which in our case is the direc- 
tion of the columns (along the z axis). In other words, 
g z (r), is the probability to find a molecule at a distance 
n ■ r from the reference molecule 

= ™ E E E to n •»■«)• ( 3 ) 

columns i— j>i 

Here N c is the number of the molecules in a column, 
is the distance between the centers of mass of molecules i 
and j, N is the total number of molecules in the system. 
/ is a binning function, i.e. / = 1 if n ■ € [r ± A/2], 
where A is the bin size (in our case the number of bins 
was 400), otherwise / = 0. V — AL x L y is the volume 
of the bin. Note that the averaging is first performed 
in each column, and then over the columns, since the 
columns can diffuse with respect to each other. 
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r (nm) r (nm) 



FIG. 6: Probability to find another molecule of the same col- 
umn at a distance r from the reference molecule. T = 400 K. 

Figurc[S]shows this distribution function for all studied 
systems at T = 400 K. It has a series of narrow peaks, 
which indicates the ordered stacking of the molecules 
in the columns. For linear side chains the peaks are 
sharp and well-defined, with a constant distance between 
molecules in the columns. It is interesting that the sep- 
aration is not affected by the length of the chain (of 
course, in the range studied). If a phenyl ring is added, 
as it is in the case of PhCi2 or branched side chains, 
the peaks broaden and shift, indicating a greater disor- 
der of the mesophase. For PhCi2 the distance between 
the columns also becomes slightly bigger (see Table IVII 
for exact numbers). Simulations at 300 K show, for all 
the systems investigated, the same distribution functions. 
This indicates that the difference of 100 K in temperature 
is not affecting the ordering properties along the columns. 
Intcrmolecular separations, which are given by the dis- 
tances between the neighboring peaks of the distribution 
function, agree well with the X-ray scattering data and 
are summarized in Table IVT1 

The second distribution function, g xyi provides infor- 
mation on the arrangement of molecules in planes per- 
pendicular to the director n, in the xy plane 

2 N N 

9xy{r) = "w^w^^ /(r ' nxnj) ' (4) 

^ ' i=0 j>i 

where V(r) = ttA(2R+5)L z is the volume of a cylindrical 
bin. 

The two-dimensional radial distribution function g xy 
is shown in Fig. [7| It has a large peak around r = 0, 
which is because of the molecules in the same column. 
The other peak is due to the in-plane nearest neighbors 
and provides us with the lattice constant of hexagonal 
(or square) arrangement of the columns. 

The spacing between in-plane (along the z-axis) near- 
est neighbors increases almost constantly with increasing 
length of the linear side chains. In the case of branched 



FIG. 7: Probability of finding the center of mass of another 
molecule at a distance r from the reference molecule (two- 
dimensional radial distribution function). The inset shows 
the large peak around r = 0. T = 400 K. 



chains, the position of the peak is almost in the same 
position as the peak belonging to the simulation of C10, 
but the peak is broader due to the disorder. For PhCi2 
side chains the peak position is between C14 and C16. 

Both the measured values and the values estimated 
from X-ray diffraction experiments are presented in Ta- 
ble and are in a fair agreement with each other. 



C. Ordering of the side chains 

To complete the description of the molecular arrange- 
ment in our systems, we also investigated the conforma- 
tional state of the side chains and the degree of their in- 
terdigitation. The conformational behavior of the alkyl 
side chains is shown in the upper panel of Fig. [SJ The 
distribution functions are given for the C12 system. The 
other systems show identical results. The distribution 
function of the dihedral CA-C2-C2-C2 (see Fig.|3]for the 
atom names), that is the dihedral connecting the core of 
the molecule with the alkyl side chains, shows that the 
gauche conformation (at around 60°) has the same prob- 
ability as the trans one (at around 180°). The distribu- 
tion functions of the dihedrals along the alkyl chains (C2- 
C2-C2-C2) and the terminating one (C2-C2-C2-C3)show 
that the trans state is more probable than the gauche 
one. For the PI1C12 system the distribution function of 
the CH-CP-CP-CH dihedral, shown in the bottom panel 
of Fig. 121 indicates that the phenyl ring tends to align 
perpendicularly to the core of the molecule. 

We have also evaluated the degree of interdigitation 
of the side chains. The results, which are shown in 
Table E] were evaluated by subtracting the in-plane 
spacing between molecules in neighboring columns from 
twice the averaged end-to-end distances of the side chains 
and a (constant) distance between the center of mass of 
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FIG. 8: Distribution function of the dihedral angles. The 
upper panel refer to C12 run and the lower one to PhCi2, 
both at T=400K. 



system 


end-to-end (nm) 


interdig. (nm) 


interdig. (%) 


C10 


1.02 


0.61 


59.8 


C12 


1.19 


0.76 


63.8 


C14 


1.35 


0.91 


67.4 


Cie 


1.51 


1.07 


70.9 


PhCia 


1.50 


1.14 


76.0 



TABLE V: Interdigitation of the side chains. 



the molecule and the atom CA (which is approximately 
0.57 nm). The longer the side chain is, the higher the 
percentage of interdigitation is. 



D. Translational dynamics 

To show that the systems have liquid- like dynamics, we 
calculated the mean square displacement of the molecular 
center of mass, along the column direction (z axis) 

MSD(t) = ([z(t)-z(0)] 2 ), (5) 

where the average is taken over all molecules. 

The resulting profiles are shown in Fig. The log-log 
profiles show that the mean-squared displacement has 1 5 
asymptotics for small t which becomes t 1 in the long- 
time limit. In fact, this type of asymptotics is expected 
for single-file diffusion |35l l86T |. when individual parti- 
cles are unable to pass each other and the particle order 
remains the same over time. In this case the motion of in- 
dividual particles requires the collective motion of many 
other particles in the same direction. In other words, in 
a one-dimensional geometry the single-particle motion is 
coupled to the density fluctuations. The fact that we can 
see the t 1 asymptotics of the mean-squared displacement 
indicates that we already see diffusion of the column as a 
whole. Note that, in an infinite one-dimensional system 



FIG. 9: log-log plot of the mean squared displacement. T = 
400 K. 



there is no diffusion, since the diffusion coefficient scales 
as 1/N. 

For T = 300 K the log-log profiles show that the mean- 
squared displacement has different asymptotics (smaller 
than t 1 / 2 ) than for T = 400 K. However, the length of 
simulation runs does not allow us to conclude wether we 
have a frozen liquid crystalline state or a solid state. 



IV. DISCUSSION AND CONCLUSIONS 

We first compare our results to the experimental data 
available in Refs. jUHE^. 

The X-ray diffraction experiments predict a columnar 
phase Dh (hexagonal arrangement of columns) at 400 K 
and D r (rectangular arrangement) at 300 K. Our simula- 
tions show that the systems at 400 K keep the hexagonal 
symmetry of the lattice during the simulation run. The 
systems at 300 K stay rectangular for short side chains; 
for C14, C16 systems the square lattice evolves into a 
hexagonal one after about 2 x 10 5 ps. 

Concerning xy spacing between the columns, the dis- 
tance increases with the length of the side chains. This 
trend seems not to be present in experimental data, 
where, for example, the spacing for C14 and C16 is iden- 
tical. The deviation between the experimental data and 
simulation results here is about 10%. 

The simulation data concerning the spacing between 
molecules in the z-direction (along the columns) are in 
fair agreement with the experimental results. Maximum 
deviation for linear chains is about 3%; for the phenyl 
case there is a deviation of about 10%. 

The order parameter obtained in simulations is sys- 
tematically higher than the one obtained experimentally. 
Of course, NMR data do not provide exact numbers for 
the order parameter, rather an estimate of it. However, 
the main problem is the extremely small size of the stud- 
ied systems, since periodic boundaries do not allow long- 
range fluctuations of the molecular order, as well as pro- 
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hibiting the presence of defects in the system. Therefore, 
calculated values of the order parameter merely empha- 
size the difference between different systems, in a quali- 
tative manner: the linear side chains provide perfect or- 
dering of the molecules in the columns; Ph— C12 systems 
show good ordering; branched side chains introduce a sig- 
nificant amount of disorder into the system, often leading 
to domains of material of different orientation. 

Finally, translational dynamics shows typical single- 
file diffusion behavior. It also indicates that one needs at 
least 10 5 ps in order to calculate correct dynamic averages 
and to equilibrate the systems. 

The main message of this work is that atomistic molec- 
ular dynamics simulations do not really provide complete 
quantitative description of the columnar mesophases be- 



cause of the very small time- and length-scales they can 
reach. Hence, their main task is to serve as an input to 
coarse-graining techniques [^3, HE E3- m order to 
devise appropriate coarse-grained representations. One 
can then treat much bigger systems on longer timescales. 
This work is in progress. 
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p,400K 
g cm" 3 


p,300K 
g cm -3 


d, 400K 
nm 


d, 300K 
nm 


h, 400K 
nm 


h, 300K 
nm 


Q, 400K 


Q, 300K 


C10-6 


1.001 


1.044 


2.54 


2.62 


0.363 


0.331 


0.85 


0.97 


CIO 


1.043 


1.087 


2.57 


2.26 


0.367 


0.363 


0.98 


0.98 


C12 


1.027 


1.074 


2.86 [23] 
2.52 [26] 
2.76 


2.38 x 6.1 °'[23] 
2.53 


0.355 [26] 
0.35 [23] 
0.366 


0.365 


0.78 [26] 
0.84 [24] 
0.98 


0.98 


C14 


1.012 


1.013 


2.66 [26J 
2.93 


2.69 


0.367 [26J 
0.367 


0.362 


0.98 


0.98 


C16 


1.002 


1.002 


2.66 [26] 
3.09 


2.87 


0.372 [26] 
0.368 


0.363 


0.98 


0.98 


Ph-C12 


1.062 


1.103 


3.4 [23] 
3.00 


2.98 [23] 
2.69 


0.35 [24] 
0.346 [23J 
0.379 


0.439 


0.93 [24] 
0.95 


0.92 



"this phase has the herringbone symmetry, see Fig. IT1 



TABLE VI: Summary of the experimental data and simulations. Results of the simulations are marked in bold. 



